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I. INTRODUCTION 



The flux of atmospheric tau neutrinos in the GeV to TeV energy range comes from both 
the intrinsic atmospheric Vj. flux and the flux due to the neutrino flavor oscillation v^-. 
The importance of understanding such a flux is twofold. First, the detection of atmospheric 
Vr flux is important for conflrming the atmospheric oscillation scenario which is 

so far established only by the z/^ disappearance measurement |lj. Second, the atmospheric 
flux is also an important background for the search of astrophysical Vr fluxes or 
exotic fluxes such as those arising from dark matter annihilations if an effective tau 
neutrino astronomy can be developed in the future. The techniques for identifying the 
tau neutrinos in the GeV to TeV energy range are discussed both in the experiments for 
atmospheric tau neutrinos and in the accelerator-based neutrino experiments 0, Q|- Due 
to growing attentions on the direct detections, it is important to investigate closely the 
flux of atmospheric tau neutrinos in such an energy range. 

The calculation of atmospheric z/^ and z/g fluxes has reached to a rather advanced stage 
The upward atmospheric z/^ flux can be calculated easily from the above flux by 
multiplying the z/^ oscillation probability arising from neutrino propagation inside the 
Earth. In such a case, the z/^ oscillations in the atmosphere is generally negligible. 

However, the atmospheric neutrino coming with a zenith angle ^ less than 90° directly arrives 
at the detector. Hence oscillations of these neutrinos in the atmosphere are precisely the 
problems one has to deal with. In this regard, an estimation of the atmospheric flux is 
given in Ref . \A for the two-flavor neutrino oscillation framework l9| . A detailed calculation 
of this flux in the same neutrino oscillation framework is given in [10|| for zenith angles 
< ^ < 60°, where the Earth curvature can be neglected in the calculation. The extension 
of such a calculation to large zenith angles is important. First of all, the path-length for the 
neutrino propagation in the atmosphere increases drastically from ^ = 60° to = 90°. As 
a result, the atmospheric flux also increases drastically in this zenith angle range by the 
enhanced z/^ oscillation probabilities. Secondly, the calculation of upward atmospheric 
Vt- flux is also important since most of the neutrino telescopes aim at detecting upward 
neutrino fluxes. In this paper, we shall extend the calculation in jl^ to the entire zenith 
angle range. 

The calculations of atmospheric flux for zenith angles 60° < < 90° and zenith 
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angles slightly larger than 90° are much more involved as we shall discuss later. Besides 
extending our previous calculation to the entire zenith angles, we also extend its validity 
from the energy range E^, > 10 GeV to the energy range E^, > 1 GeV. This improvement 
is accomplished by including the muon-decay contribution fi^ — > v^v^.e^ and /i^ z/^z/gC" 
to the intrinsic atmospheric z/^ flux, in additional to those arising from two-body vr and K 
decays. Such a contribution also generates v^- flux by z/^j — > oscillations. This part of Vr 
flux is non-negligible for E^ < 10 GeV. Furthermore, it contributes to the total z/^ flux in a 
growing percentage as the zenith angle ^ increases. 

The paper is organized as follows. In Sec. II, we introduce the method for calculating 
the intrinsic muon and tau neutrino fluxes. Particularly we outline the strategy for dealing 
with atmospheric neutrino flux for ^ > 60°. In Sec. Ill, we present the atmospheric tau 
neutrino flux taking into account the neutrino flavor oscillations. We discuss implications 
of our results in Sec. IV. 



II. THE INTRINSIC ATMOSPHERIC NEUTRINO FLUXES 
A. Intrinsic atmospheric muon neutrino flux 



We follow the approach in [ll[ for computing the flux of intrinsic atmospheric muon 
neutrinos which could oscillate into tau neutrinos. This approach computes the flux of muon 
neutrinos coming from pion and kaon decays. The method for computing muon neutrinos 
arising from muon decays will be discussed later. The flux arising from vr decays reads: 



dEdX Je Je d^E^{l - r^j Jq Xn 

X^FNn{E^,EN) X exp ( ) 0Ar(^iv), (1) 

where E is the neutrino energy, ^ is the zenith angle in the direction of incident cosmic-ray 
nucleons, r^, = m^/m^, dj, is the pion decay length in units of g/cm^, A at is the nucleon 
interaction length while A^v is the corresponding nucleon attenuation length, and (pN^Ej^) 
is the primary cosmic-ray spectrum. For the simplicity in discussions, we only consider the 
proton component of 0Ar, which is given by P] 



/ , \ -2.74 

= 1.49- fEp + 2.15- exp(-0.21v/^)j , (2) 



in units of cm~^s~^sr~^GeV~^. Since our concerned energy range for the primary cosmic 
ray flux is between 10^ GeV and 10^ GeV per nucleon (corresponding to roughly a neutrino 
energy range between 10*^ GeV and 10^ GeV), the contribution by the heavier nuclei on 
the neutrino flux is between 25% and 40% according to Fig. 7 of Ref. [8]. Hence one 
expects the eventual atmospheric z/,- flux is underestimated by 25% to 40% by considering 
only the proton component of the primary cosmic ray flux. The function Pt^{Et^,X,X') is 
the probability that a charged pion produced at the slant depth X' (g/cm^) survives to the 
depth X {> X'), F]^t,{Et,, Ejy) is the normalized inclusive cross section for A^ + air — > tc'^ + Y 
given by JjJ: 

Fr^AE., Er,) ^ EjlME!^M = _ + c_(l - xy- , (3) 



where x = E^^/Eiy, c+ = 0.92, c_ = 0.81, p+ = 4.1, and p_ = 4.8. We remark that 
c+(l — xy+ corresponds to the tt~^ production while c_(l — xY^ corresponds to the 7i~ 
production. The kaon contribution to the atmospheric z/^ flux has the same form as Eq. 
with an inclusion of the branching ratio B{K — * fiu) = 0.635 and appropriate replacements 
in kinematic factors and the normalized inclusive cross section. In particular, Fj^jxiEx, -^tv) 
can be parameterized as Eq. Q with c+ = 0.037, C-_ = 0.045, p+ = 0.87, and p_ = 3.5. 
Finally the nucleon interaction length, Aat, and the nucleon attenuation length, A^v, are both 



model dependent. A simplified approach based upon the Feynman scaling render both A 



and Ajv energy independent and Zpp = 1 — Xp/Ap = 0.263 [j^llSl, whereas a PYTHIA [l4| 
calculation give rise to an energy dependent Zpp [l^. Both results on Zpp are compared 
in Fig. [U where we have extrapolated the energy dependent Zpp{E) in Ref. [3| down to 
E = 1 GeV. The above two approaches for calculating the hadronic Z moments also give 
rise to different results for Z^^^^, Z^k, ^tvtt and Zi\ik, where the last two Z- moments are 
related to the productions of pions and kaons by the nucleon-air collisions. In this paper, 
we shall only study the Zpp dependence of the atmospheric z/^ flux (and consequently the 
atmospheric flux) since the dependencies of this flux on Z^ytt and Zjyx have been studied 
in [3]. Furthermore, compared to the Zpp case, the values of and Z^k obtained by the 
Feynman scaling do not differ signiflcantly from those obtained by the PYTHIA calculations, 
as seen from 

To proceed for calculating d^A^^ {E, ^, X) / dEdX , we note that Pt,{Et,, X, X') is given by 
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X -X'\ ( m^c dT \ 



P,(E„ X, A") = exp j . exp ^--^ _J . (4) 

where A^^ = 160 g/cw? is the pion attenuation constant, r^r is the pion hfetime at its rest 
frame, while p(T) is the atmosphere mass density at the slant depth T. For ^ < 60°, the 
curvature of the Earth can be neglected so that p(T) = Tcos^/Ziq with = 6.4 km the scale 
height for an exponential atmosphere. In this approximation, the above survival probability 
can be written as |ll| 



X-X'\ fX' 



P.(E., X, X') = exp (^-^^^ J ■ Y J ' 

where = m^c^/io/cT7r is the pion decay constant. Depending on the zenith angle, we 
apply either Eq. (0} or Eq. (0) to perform the calculations. The kaon survival probability 
Pk{Ek,X,X') has the same form as P^(i?^,X, X') except replacing with Ak and 
with ex- The two-body decay contribution to the atmospheric flux is given by the sum 
of d^X;^(E,e,X)/dEdX and d^X^(E, X)/dPdX. 

We recall that Eq. (Q) and its corresponding form in the kaon decay case only calculate 
the flux of muon neutrinos arising from two-body pion and kaon decays. To calculate the 



contribution from three-body muon decays, it is useful to first obtain the muon flux IjJ: 

'^L ' = dE^ dE^ / dX"P^(i?,X,X") 
ci-c/ Je' Je' Jo 

d^E^[l - r^j Jq Xm 

X -^FNn{E^,EN) X exp ( ) <Pn{En), (6) 

where E' and E are muon energies at slant depths X" and X respectively, while Pi_i{E, X, X") 
is the muon survival probability given by lisjl 



P^(E,X,X") =exp 



^ 1 
dT- 



(7) 



T, Jr' E{T-X'\E')p{T) 

where is the muon lifetime at its rest frame and E(T — X", E') is the muon energy at the 
slant depth T with E' the muon energy at its production point X". For the zenith angle 
^ < 60°, the above survival probability can be written as 

/X" E tf,/(E cos ^+aX cos ^) 



with = m^c^ho/ CT^ the muon decay constant and a ^ 2 MeV/g/cm^ characterizing the 
muon ionization loss in the medium [l7|. Since the muons are polarized, it is convenient to 
keep track of the right-handed and left-handed muon fluxes separately. The probability for 
a produced ^~ to be right-handed or left-handed is determined by the muon polarization 



aproduc 
HQ: 



P,{^) = ^ - 77^^^ (9) 



with X = E^/ Et, and = m^^/m^. Hence Pr,l{x) = | (1 ± Pf_i{x)) are the probabilities for 
the produced muon to be right-handed or left-handed respectively. The polarization for /x"*" 
has an opposite sign to that of fi~ . The probabilities Pr^l{x) should be inserted into Eq. © 
for obtaining four different components of the muon flux: dN'^+/dE, dN'^^/dE, dN'"+/dE, 
and dN'^^/dE. There are additional four components of the muon flux arising from the kaon 
decays. The calculation of these components proceeds in the same way as the pion decay 
case. The flux resulting from the muon flux is then given by Q| 

d'Nl^^{E,i,X) _ r .^ F,t^.SE/E,) dN^^{E,,i,X) 

dEdX ^^Je ' d,{E^,X)E^ dE, ' ^ ^ 

where dn{Efj_,X) is the muon decay length in units of g/cm^ at the slant depth X and 
F^±^^^{E/E^) is the decay distribution of fif — > u^. Precisely, in the ultra-relativistic limit, 
one has [isf 

F^-^,^{y) = go{y) + PMy)^ (H) 

with go{y) = 5/3 — Sy"^ + 4?/'^/3, gi{y) = 1/3 — 3y^ + 8y'^/3. We do not include the charm- 
hadron decay contribution to the muon neutrino flux. It is shown in Ref. jiol that charm- 
hadron decays contribute less than 5% to the overall muon neutrino flux for E^, < 10^ 
GeV. 



B. Intrinsic atmospheric tau neutrino flux 

To completely determine the atmospheric tau neutrino flux, we also need to calculate its 
intrinsic component. Since the flux of intrinsic atmospheric Ur arises from Dg decays, one 
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calculates this flux by solving the following cascade equations |l2| |: 



d4>p{E,X) ^ _tp^z ^ 
dX Xp Xp 

ay" ~ ^ ^pDsT- + ^DsDsl. 

CIA CLDs ^Ds 

dx ""'"^do: 

32 - 



(12) 



where the particle flux d(j)i{E, X) / dX denotes d Ni{E, X) / dEdX, di and Xi denote particle's 
decay and interaction length in g/cm^ respectively, and the Z moments Zij are defined by 

Z,(E,) = J^^ iE,—— _ , (13) 

with dniA-^jY{,Ei, Ej) = daiA^jviEi, Ej) / aiA^Ei). In the decay process, the scattering length 
Xi is replaced by the decay length di while driiA^jviEi, Ej)/dEj is replaced by the decay 
distribution Fi^j{Ej/Ei). In our concerned energy range, the chain of equations in ()12j) can 
be easily solved by simplifying the second equation, namely by neglecting terms (pnJ^Ds 
and Zd,Ds4>d J Xd,- One obtains 

d'N,AE,X) _ ZpoXE) {E) exp(-X/Ap)0p(E) 

d^dX l-Zpp{E) ' Ap ' ^ ^ 

We use two different values for Zpp = 1 — Xp/ Ap as shown in Fig. d To determine Zpo^ , it 
is necessary to calculate dapA^DsviEp, EaJ/dEci^. Since Ds meson is heavy enough, the 
above differential cross section is calculable using perturbative QCD j2^. In this work, the 



next-to-leading order (NLO) perturbative QCD [2l|, |22] with CTEQ6 parton distribution 
functions are employed to calculate the differential cross section of pA — > cc. To obtain 
dapA^DsviEp, EdJ /dEj:,^, we multiply the charm quark differential cross section by the 
probability factor 13% to account for the fragmentation process c ^ Ds j^]. In Fig. El we 
compare our Zpo^ to a previous result obtained by the CTEQ3 parton distribution functions 



2J|. In the latter work, the NLO pertubative QCD effects are taken into account by the K 



factor defined by 



da^^<'{E,XE)/dxE 



where da^^{E, XE)/dxE and dcr^^'-^^E, XE)/dxE are leading order and next-to- leading order 
differential cross sections for pA cc respectively, with xe = Ec/Ep. For QCD renormal- 
ization scale = rUc and the factorization scale M = 2mc, the K factor is fitted to be 



24| 

K{E,xe) = 1.36 + 0.42 In (ln(^/GeV)) 

+ (3.40 + 18.7(E/GeV)-°-^^ - 0.079 ln(E/GeV)) (16) 

We apply this K factor to our calculation with CTEQ6 parton distribution functions. Com- 
paring this result with that obtained by applying CTEQ3 parton distribution functions, one 
acquires an idea on the uncertainty of perturbative QCD approach to the charm hadron 
production cross section. It is seen from Fig. |2 that both Z moments agree well for energies 
below TeV. For E = 10 TeV, they differ by about 30%. 

Besides perturbative QCD approach, there are non-perturbative approach for comput- 
ing the charm hadron production cross section. In fact, such non-perturbative approaches 
13, are motivated to accommodate accelerator data on strange particle productions, 
which are underestimated by the perturbative QCD approach. It is desirable to apply these 
approaches to charm hadron productions. The quark-gluon-string-model (QGSM) [25;,] is a 
non-perturbative approach based upon the string fragmentation, where the model parame- 
ters are tuned to the production cross section of strange particles. The recombination-quark- 
parton-model (RQPM) |j2^ is also a phenomenological approach which takes into account 
the contribution of the intrinsic charm in the nucleon to the charm hadron production cross 
section. Detailed comparisons of these two models with perturbative QCD approach on the 



charm hadron productions are given in 27|. It is shown that perturbative QCD approach 
gives the smallest charm production Z moments. It is clear that the model dependencies on 
the charm hadron productions affect both the prompt atmospheric muon neutrino flux and 
the intrinsic atmospheric tau neutrino flux. A detailed study on the model dependencies of 
the intrinsic atmospheric tau neutrino flux is given in We shall further discuss these 
model dependencies after commenting on the Z moment Zij^,^^. 

We note that Zjy^^^ is related to the energy distributions of the Dg decays into tau neu- 
trinos. One arises from the decay Dg —>■ i^tT, and the other follows from the subsequent 
tau-lepton decay, r — >• z/^ + X. The latter contribution is calculated using the decay distri- 
butions of the decay modes r — > z/^p, r — > z/^vr, r — > z/^ai and r i't^i'i ll^. 

The uncertainty of intrinsic atmospheric Ur flux due to different approaches for Zpp is 
negligible. The main uncertainty of this flux is due to the model dependence of the Z-moment 
ZpDs- Within the perturbative QCD approach, the dependence of this flux on the parton 
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distribution functions is shown in Fig. El It is easily seen that the intrinsic atmospheric v-r 
flux is not sensitive to parton distribution functions for E < 10^ GeV. However, aX E = 10^ 
GeV, both fluxes differ by almost 50%. Incorporating the non-perturbative approaches 



for charm hadron productions [25|,|26[, the uncertainties of intrinsic atmospheric Ur flux is 
depicted in Fig. |3 It is seen that the minimal v^. flux in Ref. j3l is consistent with our z/^ 
flux calculated by perturbative QCD with CTEQ6 parton distribution functions. On the 
other hand, the maximal flux shown in Fig. |3] is almost one order of magnitude larger than 
the minimal one. This maximal flux is given by the RQPM model below 300 GeV while it 
is given by the QGSM model beyond this energy j^J]. We remark that the original minimal 
and maximal z/^ fluxes in Ref. [2^ correspond to different sets of primary cosmic ray flux, 
which is considered as one of the uncertainties for the flux. However, we have re-scaled 
these fluxes to a common cosmic ray flux, Eq. used in this paper. We also note that the 
uncertainty of intrinsic atmospheric Ur flux provided by Ref. 28[ starts aX E = 100 GeV, 
while our calculation of this flux starts at ii^ = 10 GeV. 

It is interesting to see how much the uncertainty of the intrinsic Vr flux could affect the 
determination of the flux taking into account the oscillation effect. In the next section, we 
shall study this issue with respect to the upward atmospheric flux where the oscillation 
effect is the largest. 



III. THE ATMOSPHERIC TAU NEUTRINO FLUX WITH OSCILLATIONS 



A. The Downward and Horizontal Atmospheric Tau Neutrino Fluxes 



The atmospheric tau neutrino flux can be calculated using 



dE 



^max(C) 



dX 



d^x,^(E,e,x) 



dEdX 

d'N,^{E,^,X) 



P.,^.,(i?,L(X,0) 



dEdX 



{l-P.,^uAE,L{X,0)) 



(17) 



where Pv^^v^ {E,L{X,^)) = sin^ 26'23 sin^(1.27Am3]^L/_E) is the — > Ur oscillation proba- 
bility, assuming a vanishing 6*13. We have used the notation dN^^{E,C,)/dE to denote the 
atmospheric u^- flux taking into account the oscillation effect. The unit of Am|^ is eV^ while 
L and E are in units of km and GeV respectively. A recent SK analysis of the atmospheric 
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neutrino data implies ^ 

Aml^ = (1.9 - 3.0) ■ 10"^ eV^ sin^ 2^23 > 0.9. (18) 

This is a 90% CL. range with the best fit values given by Am|]^ = 2.4 ■ 10^'^ eV^ and 
sin^ 26*23 = 1 respectively. 

1. Meson decay contributions 

Using the best fit values of neutrino oscillation parameters, we obtain atmospheric tau 
neutrino fluxes for cos^ = 0.2, 0.4, ■ ■ ■ , 1 as depicted in Fig. El This set of result is 
obtained by using an energy-independent Z moment, Zpp = 1 — Xp/Ap = 0.263 mentioned 
earlier. For the flux on the R.H.S. of Eq. (jl7j) . we only include the two-body pion and 
kaon decay contributions. The muon-decay contribution to this flux will be presented later. 
The intrinsic flux in the same equation is taken to be that calculated by perturbative 

23| . It is instructive to separately present 



QCD with CTEQ6 parton distribution functions 
the oscillated and intrinsic atmospheric Ur fluxes corresponding to the two terms on the 
R.H.S. of Eq. (jl7j) . This is done for ,^ = 0° in Fig. El We see that the oscillated and intrinsic 
atmospheric fluxes cross at -E ~ 20 GeV, indicating that the z/^ Vj. oscillation effect 
becomes important for E <2Q GeV for such a zenith angle. 

We note that the atmospheric Vr flux increases as ^ increases from 0° to 90°. There are 
two crucial factors dictating the angular dependence of such a flux. First, the atmosphere 
depth traversed by the cosmic ray particles increases as the zenith angle ^ increases. Second, 
the atmospheric muon neutrinos are on-average produced more far away from the ground 
detector for a larger zenith angle, implying a larger oscillation probability. In 

fact, the neutrino path-length dependencies on the zenith angle ^ and the neutrino energy E 
have been studied carefully by the Monte-Carlo simulation jSlj. Our semi-analytic approach 
reproduces these dependencies very well. It is found that, for = 10 GeV and cos^ = 1 
(^ = 0°), the average neutrino path-length from the z/^ production point to the ground 
detector is 14 km. The average neutrino path-length increases to 45 km and 650 km for 
£^ = 66° (cos.^ = 0.4) and C, = 90° respectively. The huge path-length of horizontal neutrinos 
makes the flux in this direction two orders of magnitude larger than the downward Ur 
flux. It is also interesting to note that the horizontal Ur flux for E approaching 1 GeV begins 
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to show oscillatory behavior. This is because, for E = 1 GeV and Am|]^ = 2.4 ■ 10~^ eV^, 
Lose = 4:E/Am'^-^ ^ 330 km which is already shorter than the average neutrino path-length 
at this zenith angle. 

We stress that our calculation procedures for cos^ > 0.5 and cos^ < 0.5 are different. In 
the former case, the curvature of the Earth can be neglected and the pion or kaon survival 
probability in the atmosphere is approximated by Eq. (0). This is the approach we adopted 
in Ref. 10|. For cos,^ < 0.5, i.e., ^ > 60°, we use Eq. (H)) for the meson survival probability. 
In this case the calculation is much more involved as the meson survival probability in 
Eq. contains an additional integration. It has been pointed out in Ref. 0| that one may 
apply Eq. ^ for calculating the path-length distribution of neutrinos for ^ > 60° so long 
as one replaces cos^ by cosefr^, where the latter is a fitted function of the former. Precisely 
speaking, by fitting the analytic calculation based upon Eq. (0) Monte-Carlo 
calculation, the relations between cos ^ and coSog ^ can be found, which are tabulated in 
Extrapolating such a relation, we find that coscg^ = 0.05 for cos^ = 0. Using this 
coSeff^ with Eq. we also calculate the atmospheric Ur flux. The result is compared with 
that obtained by the full calculation (applying Eq. (Q) as shown in Fig. [3 Both results 
agree very well. Such an agreement makes our calculation compelling and also validates the 
above extrapolation on coScg ^. 

We have so far computed the atmospheric neutrino flux with an energy independent Z 
moment, Zpp = 1 — Xp/Ap = 0.263. It is important to check the sensitivity of atmospheric 
Ur flux on the values of Zpp. We recall that different results for Zpp are shown in Fig. HI 
At energies between 10^ GeV and 10^ GeV, the values of Zpp generated by PYTHIA 
slightly depend on the energy and roughly twice larger than the value we have so far used 
for calculations. We check the effect of Zpp by calculating the atmospheric u^- flux with the 
PYTHIA-generated Zpp. The comparison of this result with the earlier one obtained by 
setting Zpp = 0.263 is shown in Fig. H for ^ = 0° and Fig. H for ^ = 90°. For ^ = 0°, two 
set of results do not exhibit noticeable difference until ii^ > 10 GeV. At E = 100 GeV, they 
differ by 45%. For ^ = 90°, two results differ by 46% at E = 1 GeV while they differ by 
29% ai E = 100 GeV. Obviously, the behavior of Zpp is one of the major uncertainties for 
determining the atmospheric u^- flux. 
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2. Muon-Decay contributions 



We have stated that the muon-decay contributions to i/^ is non-neghgible for neutrino 
energies less than 10 GeV. Such i^^'s can oscillate into z/^'s during their propagations in the 
atmosphere. The calculation of such a flux according to Eqs. (0) and (jlOp is rather involved. 
However, a simple approximation as presented below gives a rather accurate result for this 
flux. 

The calculation of z/^ spectrum due to muon decays requires the knowledge of muon 
polarizations. The muon polarization however depends on the ratio of muon momentum 
to the momentum of parent pion or kaon as indicated by Eq. Q. It is straightforward to 
calculate the average muon polarization at any slant depth X provided the energy spectrum 
of the parent pion or kaon is known at that point. For the downward case = 0°) , it is known 
from the previous section that the muons are most likely produced at around 14 km from the 
ground detector. At that point, the pion and kaon fluxes can be approximately parameterized 
as MEn) = 10-^-^^ ■ E-^-°^ and (pKiEx) = W'^'^^ ■ E^^'^^ in units of cm-^s-^sr-^GeV-^ 
for meson energies between 1 and few tens of GeV. We do not distinguish 7r~(i^'~) from 
7r+(iC+) in the above fittings. Although the spectra are charge dependent, the resulting 
absolute values of /i^ polarization and polarization differ by only 10% for up to few 
tens of GeV flj|. From Eq. 0, and the above pion and kaon spectra, we obtain (-P^-) = 0.35, 
(P^) = 0.95. Therefore /i~ coming from the tt~ decays are 67% right-handed polarized and 
33% left-handed polarized. On the other hand, fi~ coming from K~ decays are 98% right- 
handed polarized and only 2% left-handed polarized. The muons produced by meson decays 
lose energies before they decay into neutrinos. The decay distribution for /i~ —>■ is given 
by Eq. ()11|) . The average momentum fraction (y) of muon neutrinos are 0.3 and 0.4 from 
decays of right-handed /i~ and left-handed fi~ respectively. Following a similar procedure, 
one can determine the polarization and decay distributions of fi^. Finally, to calculate 
the spectrum of muon neutrinos arising from muon decays, we use the approximation of 
replacing F^±^^^{E/E^) with 6{E/E^ - {y)) in Eq. 

To check the validity of the above approximation, we compare our result on the fraction 
of muon decay contribution to the overall flux with that given by Ref. yi for cos^ = 0.4, 
i.e., ^ = 66°. At this zenith angle, most of the muons are produced roughly 45 km from the 
detector. The pion and kaon fluxes at this point are fitted to be ^^^(-E'Tr) = 10~^'^^ ■ E~^'^^ and 
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^j^{Ek) = 10-^-^^ ■ E]^'^-^^ in units of cm-^s-^sr-^GeV-^ This gives rise to (P^.) = 0.34, 
{P^-) = 0.94. Following the procedure in the downward case, we obtain the muon neutrino 
flux from the muon decays. At E = 1 GeV, the fraction of muon-decay contributions to the 
overall flux is 44% while the fraction decreases to 17% at P = 10 GeV. In Ref. Q], the 
corresponding fractions are 47% and 18% respectively. Both set of fractions agree rather 
well. 

Since our approximation works well for calculating the muon-decay contributions to the 
atmospheric z/^ flux, we proceed to calculate the resulting atmospheric Ur flux with Eq. (fTTj) . 
Speciflcally we only need to include the flrst term on the R.H.S. of Eq. (fT7|) because the 
second term has already been included in the two-body decay contribution. In Fig. 
those atmospheric Ut fluxes resulting from oscillations of i/^'s generated by both two- and 
three-body decays (muon decays) are compared with those resulting from the oscillations of 
z/^'s generated only by two-body decays. As expected, the three-body decay contribution 
is non-negligible for E < 10 GeV. Quantitatively, for ^ = 0° and E = 1 GeV, 24% of the 
total atmospheric flux is from the oscillations of i/^'s originated from the muon decays. 
At = 10 GeV, only 2.9% of the total atmospheric Uj. flux comes from the same source. 
For ^ = 66°, the three-body decay contribution gives rise to 36% and 8.9% of the total 
atmospheric z/^ flux at E = 1 GeV and E = 10 GeV respectively. Finally, for C, = 90°, 
the three-body decay contribution to the total atmospheric z/^ flux is most signiflcant. It 
contributes to 53%, 46%, and 39% of the total atmospheric flux at _E = 1 GeV, 10 
GeV and 20 GeV respectively. To calculate the three-body decay contribution to z/^ flux at 
C, = 90°, we have used Eq. (0) for the meson survival probability with coSeg^ = 0.05 and a 
overall factor C ~ 1.40 to flx the normalization of the flux 



B. The Upward Atmospheric Tau Neutrino Flux 

The upward atmospheric Ur fluxes are enhanced compared to those of other directions 
since the average neutrino path lengths in such case are larger. Therefore the observations 
of astrophysical tau neutrinos in upward directions are subject to more serious background 
problems. However, observing the atmospheric tau neutrinos is interesting in its own right. 
The atmospheric tau neutrino flux for cos^ = —0.2 is shown in Fig. ^2 The effect of 
oscillation is evident for below TeV energies. This is seen from the crossing point 
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of intrinsic and oscillated atmospheric Vr fluxes at ~ 700 GeV. The atmospheric Vr flux 
shows oscillatory behavior for E < IQ GeV. For cos.^ < —0.2, such an oscillatory behavior 
is even more signiflcant. In such a case, it is more practical to study the averaged flux. We 
average the atmospheric v^- flux for the zenith angle range —1 < cos.^ < —0.4, as shown in 
Fig. ^1 Due to uncertainties of the intrinsic atmospheric flux as discussed in Ref . |28| , 
the atmospheric Vr flux taking into account the oscillation effect also contains uncertainties 
beginning at a few hundred GeV's. In the same flgure, we also plot the corresponding 
atmospheric flux. The and Vr fluxes are comparable for < 40 GeV. In such a 
case, the footprint of Vr might be identified by studying the energy spectra of shower events 
induced by neutrino interactions P|. hi E = 10^ GeV, the z/^ flux is approximately 30 times 
larger than the maximal flux. We note that the maximal and minimal Vr fluxes begin to 
differ a,t E = 500 GeV. Ai E = 1 TeV, the maximal flux is 3 times larger than the minimal 
one. The ratio of maximal flux to the minimal one increases to 14 at = 10 TeV. We 
remark that the upward atmospheric Ur flux is also calculated in Ref. with sin^ 2^23 = 1 
and Am|^ = 10~^, 10~^'^, 10"'^ eV^ respectively. Here we have done the calculation with the 
best flt value of sin^ 26'23 and Amg^ taken from |])| . Furthermore we include the contribution 
of intrinsic atmospheric Vr flux and its associated uncertainties. 



IV. DISCUSSION AND CONCLUSION 



The understanding of atmospheric flux is important for exploring the tau neutrino 
astronomy [3J3- As mentioned earlier, an estimation of the atmospheric Vj- flux has been 
given in Ref. [2] while a detailed calculation of this flux for zenith angles < ^ < 60° is given 
in Ref. In these works, comparisons of the galactic-plane flux with the atmospheric 
v-r flux are also made for illustrating the possibility of the tau neutrino astronomy. Now that 
we have obtained a complete result of the atmospheric Vr flux for the entire zenith angle 
range, we compare this flux with two astrophysical fluxes: the galactic-plane tau neutrino 
lux just mentioned and the cosmological flux due to neutralino dark matter annihilations 
3^ . The comparison is depicted in Fig. 1131 where the flux of galactic-plane tau neutrinos is 
taken from the calculation of Ref. jlfll. One can see that the galactic-plane Vr flux dominates 
over the downward = 0°) atmospheric v^- flux for E greater than a few GeV. Hence, in this 
direction, it is possible to observe the flux of galactic-plane tau neutrinos in the GeV energy 
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range. For near horizontal directions, the atmospheric Vr flux grows rapidly with zenith 
angles. Therefore, for = 90°, the energy threshold for galactic-plane tau neutrino flux to 
dominate over its atmospheric counterpart is pushed up to ~ 100 GeV. We further see 
that the galactic-plane z/,- flux does not dominate the upward atmospheric z/^ background 
(— 1 < cos^ < —0.4) until E = 500 GeV. However, it is noteworthy that, in the muon 
neutrino case, galact ic-p lane neutrino flux is overwhelmed by the atmospheric background 
until Ely > 10^ GeV |33|. Such a difference between z/^ and shows the promise of the tau 
neutrino astronomy in the GeV energy range as pointed out in . From Fig. it is 

also clear that the atmospheric Ur flux is a non-negligible background to the cosmological 
tau neutrino flux due to neutralino dark matter annihilations [32]. In fact, two fluxes are 
comparable in the downward direction while the atmospheric flux dominates in horizontal 
and upward directions. 

In summary, we have presented a semi-analytical calculation on the atmospheric flux 
in the GeV to TeV energy range for downward, upward, and horizontal directions. The 
atmospheric z/^ flux at ,^ = 90° is two orders of magnitude larger than the corresponding flux 
at ^ = 0° for 1 < E/GeV < 10. On the other hand, the fluxes with zenith angles between 
and 90 degrees merge for E > 700 GeV, provided that the intrinsic atmospheric Uj. flux 
is calculated with perturbative QCD. Should one adopt a non-perturbative model for the 
intrinsic z/^ flux, the resulting z/,- fluxes on Earth at different zenith angles would merge at an 
energy lower than 700 GeV. We have observed that the upward atmospheric z/,- fluxes show 
oscillatory behaviors. For the averaged flux with —1 < cos^ < —0.4, the atmospheric I'r 
flux is found to be comparable to the atmospheric z/^ flux for E < AO GeV. The comparison 
of this flux with the horizontal atmospheric z/^ flux is also interesting. Two fluxes are in 
fact comparable for E < 10 GeV. This shows that the z/^ oscillation is already quite 

signiflcant in the horizontal direction for such an energy range. Nevertheless, the upward 
atmospheric Vr flux takes over from E > 10 GeV until E' ~ 3 TeV where two fluxes merge 
again. Concerning the uncertainties in our calculations, we have studied the dependencies 
of atmospheric Vr flux on the Z moment Zpp for representative zenith angles ^ = 0° and 
^ = 90°. We have also discussed in detail the uncertainty of intrinsic atmospheric z/^ flux due 
to different models for charm hadron productions. The consequence of such a uncertainty 
on the determination of oscillated Vj. flux is studied as well. Concerning the technique for 
calculating the atmospheric Vr flux from large zenith angles, we have verifled the validity 
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of using coseff^ in Eq. (0) to calculate the atmospheric Vr flux for ^ > 60°. In particular, 



we have extrapolated the results in Ref. 



to = 90° and demonstrate that the choice 



COSeff(e = 90°) 

using Eq. (jH). 



0.05 reproduces well the atmospheric Vr flux obtained by a full calculation 



Acknowledgements 

This work is supported by the National Science Council of Taiwan under the grant number 
NSC 93-2112-M-009-001. 



[1] For a recent result, see Y. Ashie et al. [Super-Kamiokande Collaboration], Phys. Rev. Lett. 
93, 101801 (2004) |arXiv :hep-ex/04 04034l . 



[2] 
[3] 
K 

[5] 
[6] 

[7] 
[8] 

[9 

[10 

[11 
[12 
[13 
[14 
[15 



H. Athar, Mod. Phys. Lett. A 19, 1171 (2004). 

H. Athar and C. S. Kim, Phys. Lett. B 598, 1 (2004) |arXiv:hep-ph/0407182| . 

For a recent review, see G. Bertone, D. Hooper and J. Silk, Phys. Rept. 405, 279 (2005) 

|arXiv:hep-ph/0404175| . 

T. Stanev, Phys. Rev. Lett. 83, 5427 (1999) |arXiv:astro -ph/9907018|. 



A. Zalewska [ICARUS Collaboration], Acta Phys. Polon. B 35, 1949 (2004). 
J. Campagne, Eur. Phys. J. C 33, S837 (2004). 

T. K. Gaisser and M. Honda, Ann. Rev. Nucl. Part. Sci. 52, 153 (2002) 
|arXiv:hep-ph/0203272| . 

For an estimation of the same flux in the three-flavor neutrino oscillation framework, see 
H. Athar, |arXiv:hep-ph/0411303| 

H. Athar, F. F. Lee and G. L. Lin, Phys. Rev. D 71, 103008 (2005) |arXiv:hep-ph/0407183| . 
T. K. Gaisser, Astropart. Phys. 16, 285 (2002) |arXiv:astro-ph/0104327| . 

Cosmic Rays and Particle Physics, T. K. Gaisser, Cambridge University Press (1992). 

P. Lipari, Astropart. Phys. 1, 195 (1993). 

T. Sjostrand, Comput. Phys. Commun. 82, 74 (1994). 

M. Thunman, G. Ingelman and P. Gondolo, Astropart. Phys. 5, 309 (1996) 



|arXiv:hep-ph/9505417| . 



16 



[16] V. Agrawal, T. K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D 53, 1314 (1996) 
|arXiv:hep-ph /95094'23l . 
1] B. Rossi, High Energy Particles (Prentice Hall, Englewood Cliffs, NJ, USA, 1952). 
3] S. M. Barr, T. K. Gaisser, P. Lipari and S. Tilav, Phys. Lett. B 214, 147 (1988). 
3] G. Barr, T. K. Gaisser and T. Stanev, Phys. Rev. D 39, 3532 (1989). 
3] L. Pasquah and M. H. Reno, Phys. Rev. D 59, 093003 (1999) |arXiv:hep-ph /9811268|. 
1] P. Nason, S. Dawson and R. K. Ellis, Nucl. Phys. B 327, 49 (1989) [Erratum-ibid. B 335, 260 
(1990)]. 

2] M. L. Mangano, P. Nason and G. Ridolfi, Nucl. Phys. B 373, 295 (1992). 
3] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. Nadolsky and W. K. Tung, JHEP 0207, 

012 (2002) |arXiv:hep-ph/0Mri95| . 
4] L. Pasquah, M. H. Reno and I. Sarcevic, Phys. Rev. D 59, 034020 (1999) 
[arXiv:hep-ph/9806428, . 

5] A. B. Kaidalov and O. I. Piskunova, Z. Phys. C 30, 145 (1986); L. V. Volkova, W. Fulgione, 

P. Galeotti and O. Saavedra, Nuovo Cim. C 10, 465 (1987). 
3] E. V. Bugaev, A. Misaki, V. A. Naumov, T. S. Sinegovskaya, S. I. Sinegovsky and N. Taka- 

hashi, Phys. Rev. D 58, 054001 (1998) arXi v:hep-p h/9803488|. 
7] C. G. S. Costa, Astropart. Phys. 16, 193 (2001) |arXiv:hep-ph/0010306| . 
3] C. G. S. Costa, F. Halzen and C. Salles, Phys. Rev. D 66, 113002 (2002) 

|arXiv:hep-p h70104039 . 
3] B. A. Li, Phys. Rev. D 52, 5165 (1995) >rXiv:hep-ph/9504304| . 
3] C. G. S. Costa and C. Salles, |arXiv:hep-ph/01Q527T| 

1] T. K. Gaisser and T. Stanev, Phys. Rev. D 57, 1977 (1998) |arXiv:astro-ph/9708146| . 
I] D. Elsaesser and K. Mannheim, Astropart. Phys. 22, 65 (2004) 'ar Xiv:astro-ph/0405347| . 
3] H. Athar, K. Cheung, G. L. Lin and J. J. Tseng, Eur. Phys. J. C 33, S959 (2004) 
|arXiv:astro-ph/0311586| . 



17 



0.7 



0.6 - 



0.5 - 



UJ 
N 



0.4 - 



^ 0.3 h 



-1 1 1 r 



Our fitting 
Thunman et al. 

Z =0.263 

pp 




0.2 - 



0.1 - 



0.0 



2 3 

Log^ JE / GeV] 



FIG. 1: The comparison of Zpp obtained by assuming the Feynman scahng and that obtained 
by PYTHIA 'l^. Our extrapolation of the latter result is also shown. 
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FIG. 2: The Z moment Z^b, obtained by perturbative QCD with CTEQ3 and CTEQ6 parton 
distribution functions respectively. 
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FIG. 3: The comparison of intrinsic atmospheric Vr fluxes calculated by perturbative QCD with 
CTEQ3 and CTEQ6 parton distribution functions respectively. 
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FIG. 4: The model dependencies of intrinsic atmospheric flux. The minimal flux from Ref. 2^ 
is given by perturbative QCD. The maximal flux from the same reference is given by RQPM model 
for E < 300 GeV, and by QGSM model for E > 300 GeV. 
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E (GeV) 

FIG. 5: The atmospheric Vr Aux for cos^ = 0, 0.2, 0.4, 0.6, 0.8 and 1 (from top to bottom) with 
sin^ 26123 = 1, Am|i = 2.4 • lO'^ eV^ and Zpp = 0.263. 
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FIG. 6: The downward atmospheric Ur flux (solid Hne) as the sum of its oscillated (dotted line) 

and intrinsic (dashed line) components. 



23 




10° 10' 10' 10' 10^ 



E (GeV) 

FIG. 7: The comparison of atmospheric v-,- flux obtained by using coSeff ^ and that obtained by the 
full calculation for ^ = 90°. 



24 



Z =0.263 

pp 

Z (E) 

pp^ ' 




10 



■13 



^ = 01 ( Am\^= 2.4 . 10"'eV', sin'2e„ = 



23 



10' 



10^ 



10^ 
E (GeV) 



lO"- 



10^ 



FIG. 8: The comparison of atmospheric Vr fluxes calculated from a constant Z, 



energy-dependent Z^r, [15|] for ^ = 0°. 
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FIG. 9: The comparison of atmospheric Vr fluxes calculated from a constant Z, 



energy-dependent Z^r, [15|] for ^ = 90°. 
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FIG. 10: The comparisons of atmospheric Vr fluxes resulting from the oscillations of i/^'s generated 
by two-body and three-body decays with those resulting from the oscillations of f^'s generated 
only by two-body decays. The comparisons are made for three zenith angles, cos ^ = 0, 0.4, and 1 
(from top to bottom). 
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FIG. 11: The atmospheric Uj- flux (sohd line) as a sum of its oscillated (dotted line) and intrinsic 
(dashed line) components for cos^ = —0.2 with sin^ 2^23 = 1, ^^nii = 2.4 • 10~^ eV^, and Zpp = 
0.263. 
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FIG. 12: The atmospheric Uj- flux averaged for — 1 < cos^ < —0.4. The uncertainty of this flux 
due to the uncertainty of intrinsic atmospheric v-r flux is shown. The corresponding flux is also 
plotted, which is equal to the Vr flux for 1 < E/GeV < 10. 
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FIG. 13: The comparison of atmospheric Uj- fluxes with the galactic-plane tau neutrino flux 
and the tau neutrino flux due to the neutralino dark matter annihilations 3^. We have included 
downward (cos^ = 1), horizontal (cos^ = 0) and upward (—1 < cos^ < —0.4) atmospheric Vr 
fluxes for the comparison. At = 1 GeV, the horizontal flux is the largest whereas the downward 
flux is the smallest. 
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